classdef PeakFitting < DataProcess.Fitting.AbstractFitting
    %PEAKFITTING Summary of this class goes here
    %   Detailed explanation goes here
    
    properties
    end
    
    methods
        function obj = PeakFitting(data_x, data_y)
            %PEAKFITTING Construct an instance of this class
            %   Detailed explanation goes here
            obj@DataProcess.Fitting.AbstractFitting('PeakFitting', data_x, data_y);
        end
        
        function res = LorentzianFitting(obj, varargin)
            p=inputParser;
            p.addParameter('k', 1, @(x) x>0);
            p.addParameter('func', 'direct', @(x) ismember(lower(x), {'direct', 'sqrt'}));
            p.addParameter('width_factor', 2.0, @(x) x>0.0);
            p.addParameter('background', 0.0, @isnumeric);
            p.parse(varargin{:});
            
            idx = p.Results.k; 
            data_x = obj.data_x; data_y = obj.data_y;
            [d.peak, d.loc, d.width, d.pro] = findpeaks(data_y, data_x, 'Npeak', idx, 'SortStr', 'descend');
            pk_k = d.peak(idx); width_k = d.width(idx); loc_k = d.loc(idx);
            res.direct_info = d; 
            
            [xData, yData] = prepareCurveData(data_x, data_y);
            
            q = p.Results.width_factor;
            obj.exclude_index = find(data_x<d.loc-q*d.width | data_x>d.loc+q*d.width);
            obj.keep_index    = find(data_x>=d.loc-q*d.width & data_x<=d.loc+q*d.width);
            excludedPoints = excludedata( xData, yData, 'Indices', obj.exclude_index);

            b = p.Results.background;
            switch lower(p.Results.func)
                case 'sqrt'
                    lineshapeStr = 'a*g./sqrt(((x-x0).^2+g^2)) + b';
                case 'direct'
                    lineshapeStr = 'a*g^2/((x-x0).^2+g^2) + b';
            end
            % Set up fittype and options.
            ft = fittype(lineshapeStr, 'independent', 'x', 'dependent', 'y' );            
            opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
            opts.Display = 'Off';
            opts.Lower =      [0.5*(pk_k-b)  0.5*b  0.5*width_k  0.5*loc_k];
            opts.StartPoint = [1.0*(pk_k-b)  1.0*b  0.5*width_k  1.0*loc_k];
            opts.Upper =      [1.5*(pk_k-b)  1.5*b  1.5*width_k  1.5*loc_k];
            opts.Exclude = excludedPoints;

            % Fit model to data.
            [fitting.result, fitting.gof] = fit( xData, yData, ft, opts );
            res.fitting = fitting;
            res.center = fitting.result.x0;
            res.width = fitting.result.g;
            res.amplitude = fitting.result.a;
            res.background = fitting.result.b;
            res.rsquare = fitting.gof.rsquare;
            res.function = fitting.result;
            
            obj.result = res;
        end
        
        function plot_fitting_curve(obj, varargin)
            p=inputParser;
            p.addOptional('xList', obj.data_x, @isnumeric);
            p.addParameter('handle', [], @(x) isa(x, 'matlab.graphics.axis.Axes') || isempty(x));
            p.parse(varargin{:});
            
            if isempty(p.Results.handle)
                ax = axes;
            else
                ax = p.Results.handle;
            end
            plot(ax, p.Results.xList, obj.result.function(p.Results.xList), 'r-', ...
                     obj.data_x(obj.keep_index), obj.data_y(obj.keep_index), 'ro', ...
                     obj.data_x(obj.exclude_index), obj.data_y(obj.exclude_index), 'bx');
        end
    end
end

